function ret = ldp(img, mode)
% Input:
%   img - input image M*N
%   mode -    ''  -> output M*N*32 (values in 0 or 1 in each direction)
%             's' -> output M*N*4  (values in 0..255 in each direction)
%   (default) 'sh'-> output 256*4  (bins stand for 0..255 in each direction)
% Output:
%   ldp feature (2rd) in each pixel.
%
% Description & Example:
% 1) Weight:
%  A B C    1  8  32
%  D E F -> 2     64
%  G H I    4 16 128
% 2) Directions:
% >>  a = magic(5);
%   a = 
%   17  24  1   8   15
%   23  5   7   14  16
%   4   6   13  20  22
%   10  12  19  21  3
%   11  18  25  2   9
% At point (3,3):
% DIRECTION 1 : xor(sgn(A-B),sgn(C-D)) for all 8 C
%   17  24  1   8   15
%   D   C   7   14  16      1  0  0
%   4   B   A   20  22      0     0
%   10  12  19  21  3       0  0  0
%   11  18  25  2   9
% DIRECTION 2 : xor(sgn(A-B),sgn(C-D)) for all 8 C
%   17  D   1   8   15
%   23  C   B   14  16      1  0  0
%   4   6   A   20  22      0     0
%   10  12  19  21  3       0  0  0
%   11  18  25  2   9
% DIRECTION 3 : xor(sgn(A-B),sgn(C-D)) for all 8 C
%   D   24  1   8   15
%   23  B,C 7   14  16      1  1  0
%   4   6   A   20  22      1     0
%   10  12  19  21  3       0  0  0
%   11  18  25  2   9
% DIRECTION 4 : xor(sgn(A-B),sgn(C-D)) for all 8 C
%   17  24  D   8   15
%   23  C   7   B   16      1  0  0
%   4   6   A   20  22      0     1
%   10  12  19  21  3       0  0  0
%   11  18  25  2   9
%
% the result for a(3,3) is 1; 1; 11; 65;
%
if nargin < 2
    mode = 'sh';
end
% S - 32 bit per pixel
img = double(img);
[P, Q] = size(img);
fx = uint8(img(2:end-1, 2:end-1) - img(2:end-1, 1:end-2) > 0); % 0
fy = uint8(img(2:end-1, 2:end-1) - img(1:end-2, 2:end-1) > 0); % 90
fd = uint8(img(2:end-1, 2:end-1) - img(1:end-2, 1:end-2) > 0); % 45
fu = uint8(img(2:end-1, 2:end-1) - img(1:end-2, 3:end) > 0); % 135
s = zeros(P-4, Q-4, 32);

dx = [-1 -1 -1  0 0  1 1 1];
dy = [-1  0  1 -1 1 -1 0 1];
% split image into rectangles?
% bx = [ 3  2  3  3 3  3 2 3 ...
%        2  1  2  2 2  2 1 2 ...
%        3  2  3  3 3  3 2 3 ...
%        3  2  3  3 3  3 2 3
%     ];
% by = [ 2  2  2  1 1  2 2 2 ...
%        3  3  3  2 2  3 3 3 ...
%        3  3  3  2 2  3 3 3 ...
%        3  3  3  2 2  3 3 3
%     ];
for i = 1:8
    s(:, :, i) = bitxor(fx(2:end-1, 2:end-1), fx(2+dy(i):end-1+dy(i), 2+dx(i):end-1+dx(i)));
    s(:, :, i+8) = bitxor(fy(2:end-1, 2:end-1), fy(2+dy(i):end-1+dy(i), 2+dx(i):end-1+dx(i)));
    s(:, :, i+16) = bitxor(fd(2:end-1, 2:end-1), fd(2+dy(i):end-1+dy(i), 2+dx(i):end-1+dx(i)));
    s(:, :, i+24) = bitxor(fu(2:end-1, 2:end-1), fu(2+dy(i):end-1+dy(i), 2+dx(i):end-1+dx(i)));
end

if strcmpi(mode, 'sh') || strcmpi(mode, 's')
    v1 = bin2dec(s(:,:,1:8));
    v2 = bin2dec(s(:,:,9:16));
    v3 = bin2dec(s(:,:,17:24));
    v4 = bin2dec(s(:,:,25:32));
    if strcmpi(mode,'sh')
        % create histogram
        h1 = hist(v1(:),[0:255]);
        h2 = hist(v2(:),[0:255]);
        h3 = hist(v3(:),[0:255]);
        h4 = hist(v4(:),[0:255]);
        ret = [h1 h2 h3 h4];
    else
        ret(:,:,1) = v1;
        ret(:,:,2) = v2;
        ret(:,:,3) = v3;
        ret(:,:,4) = v4;
    end
else
    ret = s;
end

end

function val = bin2dec(data)
[M,N,P] = size(data);
val = zeros(M,N);
k = 1;
for i = 1:P
    val = val + data(:,:,i)*k;
    k = k*2;
end
end
